Magnetic Field Diagnostics Based on Far-Infrared Polarimetry: 
Tests Using Numerical Simulations 



Fabian Heitsch^ 
Ellen G. ZweibeP 
Mordecai-Mark Mac Low'^ 
Pakshing Li^ 
Michael L. Norman'* 

ABSTRACT 

The dynamical state of star-forming molecular clouds cannot be understood without deter- 
mining the structure and strength of their magnetic fields. Measurements of polarized far-infrared 
radiation from thermally aligned dust grains are used to map the orientation of the field and es- 
timate its strength, but the accuracy of the results has remained in doubt. In order to assess the 
reliability of this method, we apply it to simulated far-infrared polarization maps derived from 
three-dimensional simulations of supersonic magnetohydrodynamical turbulence, and compare 
the estimated values to the known magnetic field strengths in the simulations. We investigate the 
effects of limited telescope resolution and self-gravity on the structure of the maps. Limited ob- 
servational resolution affects the field structure such that small scale variations can be completely 
suppressed, thus giving the impression of a very homogeneous field. The Chandrasekhar-Fcrmi 
method of estimating the mean magnetic field in a turbulent medium is tested, and we suggest 
an extension to measure the rms field. Both methods yield results within a factor of 2 for field 
strengths typical of molecular clouds, with the modified version returning more reliable estimates 
for slightly weaker fields. However, neither method alone works well for very weak fields, missing 
them by a factor of up to 150. Taking the geometric mean of both methods estimates even the 
weakest fields accurately within a factor of 2.5. Limited telescope resolution leads to a system- 
atic overestimation of the field strengths for all methods. We discuss the effects responsible for 
this overestimation and show how to extract information on the underlying (turbulent) power 
spectrum. 

Subject headings: ISMiClouds, Turbulence — ISM:Kinematics and Dynamics — ISM:Magnetic Fields — 
polarization 
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1. Motivation 

The significance of magnetic fields in the dy- 
namics of molecular clouds, and in star formation 
itself, is still uncertain. 

Theory and modelling have demonstrated that 
the strength of magnetic fields can have a pro- 
found influence on the processes which lead to 
the formation of stars. There is a critical ratio 
of magnetic flux to mass above which magnetic 
fields can prevent gravitational collapse (magnet- 
ically subcritical) and below which they cannot 
(magnetically supercritical; e.g. Mouschovias & 
Spitzer 1976). Redistribution of magnetic flux by 
ambipolar diffusion in a magnetically subcritical 
region can lead to loss of support against self- 
gravitation, and eventually to low mass star for- 
mation (see Shu, Adams, & Lizano 1987, Myers & 
Goodman (1988), Porro & Silvcstro (1993), Ciolek 
& Mouschovias 1993, 1994, 1995, Safier, McKee, 
& Stahler 1997, Ciolek & Konigl 1998). How- 
ever, magnetically siibcritical regions do not ap- 
pear to reproduce observations of molecular cloud 
cores (Nakano 1998), nor can this scenario ex- 
plain the large fraction of cloud cores contain- 
ing protostellar objects (Ward-Thompson, Motte, 
& Andre 1999). It is possible that magnetically 
supercritical regions can be supported by turbu- 
lence (Bonazzola et al. 1987, McKee & Zweibel 
1995), which, to support a highly supercritical re- 
gion, must be quite nonlinear (Myers & Zweibel 
2001). Recent 3D simulations (Heitsch, Mac Low, 
& Klessen 2001, hereafter HMK, Ostriker, Stone, 
& Gammie 2001) demonstrate that strong turbu- 
lence can provide large-scale support against col- 
lapse, though it cannot prevent collapse on small 
scales. The degree to which clouds are supported 
at all, of course, depends on their lifetimes and 
on the star formation rate (Ballesteros-Paredes, 
Hartmann, & Vazquez-Semadeni 1999, Elmegreen 
2000). 

Magnetic fields affect other dynamical pro- 
cesses in clouds besides gravitational collapse. 
In a strong magnetic field, weakly compressible 
turbulence is anisotropic (Sridhar & Goldreich 
1994; Goldreich & Sridhar 1995, 1997) and en- 
ergy dissipation is relatively more in waves and 
less in shocks (Smith, Mac Low, & Zuev 2000), 
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although the overall rate of energy dissipation is 
not strongly dependent on the fieldstrength (Mac 
Low 1999). Magnetic fields may also play a role in 
coUimating molecular outflows, and in transferring 
their momentum to the ambient medium. 

The actual strength of magnetic fields in molec- 
ular clouds will ultimately determine which theo- 
retical picture is correct, so observations of field 
strengths are crucial. The integrated line of sight 
component of the field, weighted by the density 
of the tracer species, can be measured through 
the Zeeman effect. However, Zeeman mapping is 
time constuning and rcriuircs high sensitivity and 
the presence of particular tracers (most commonly 
OH). Moreover, this technique does not probe the 
field in the plane of the sky. For all of these 
reasons, it is desirable to have a complementary 
method of mapping the field. 

Following the first detection by Cudlip et al. 

(1982), it has been shown that the polarization of 
the far infrared thermal radiation emitted by mag- 
netically aligned dust grains can be used to map 
the orientation of the magnetic field on the plane 
of the sky (see Hildebrand et al 2000 for methodol- 
ogy). This has been done, at various wavelengths, 
for a number of nearby clouds (Gonatas et al. 
1990, Jarrett et al 1994, Dotson 1996, Rao et al. 

1998, Schleuning 1998, Glenn, Walker, & Young 

1999, Dotson et al. 2000, Schleuning et al. 2000, 
Vallee, Bastien, & Greaves 2000, Ward-Thompson 
et al 2000) and near the Galactic Center (Werner 
et al. 1988, Hildebrand et al. 1990, 1993, Dowell 
1997, Novak et al. 1997,2000). Far infrared polar- 
ization maps display the morphology of the field 
relative to other structures in the cloud, and can 
also be used to estimate the strength of the mean 
field according to a dynamical method originally 
proposed by Chandrasckhar & Fermi (1953, here- 
after CF). Applications of the CF- method gener- 
ally suggest fieldstrengths in the milligauss range 
or above. Such fieldstrengths are larger than what 
is typically measured by the Zeeman effect (e.g. 
Glenn et al. 1999, Lai et al. 2000). 

The CF-mcthod is subject to errors arising from 
line of sight and angular averaging, and rests on 
the assumptions of equipartition between turbu- 
lent kinetic and magnetic energy, and isotropy of 
fluid motions. Spatial averaging, by smoothing 
the maps, can give misleading impressions about 
the magnetic field morphology. However, since the 
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information carried in the maps is unique, it is dif- 
ficult to test the magnitude of these effects with 
astronomical observations. 

Numerical simulations of turbulent, magnetized 
molecular clouds offer the means to calibrate the 
accuracy of polarization maps and develop new 
techniques to analyze them. The actual strength 
and structure of the field are known at all grid- 
points and at selected times, as are the gas veloc- 
ity and density. It is possible to create synthetic 
polarization maps, analyze them as though they 
were astronomical data, and check the accuracy 
of the results. That is the subject of this paper. 
Ostriker, Stone, & Gammie (2001) have used a set 
of simulations to calibrate the CF-method in this 
manner. 

In section 2, we describe the numerical models 
and the method by which we generate polariza- 
tion maps. Most of the results relevant to polar- 
ization maps arc in section 3, in which we discuss 
the morphology shown in the maps, implement the 
CF-method, devise an alternative to it, and show 
what can be learned about the spectrum of mag- 
netic field fluctuations. Section 4 is a summary 
and discussion. 

2. Models and Methods 

2.1. Models of 3D-MHD-turbulence 

We base our investigation on full 3D models 
of driven MHD-turbulence in a cube with peri- 
odic boimdary conditions, simulating a portion of 
the interior of a molecular cloud. We chose an 
isothermal equation of state, because the cooling 
times arc much shorter than the dynamical times 
at the high densities typical of molecular clouds. 
We performed the simulations at 128^ and 256^ 
grid zones resolution using ZEUS-3D, a well-tested 
Eulerian finite-difference code (Stone & Norman 
1992a, 1992b, Clarke 1994) with second-order ad- 
vcction and a von Neumann artificial viscosity to 
capture shocks. The MHD induction equation 
is followed using the method of consistent trans- 
port along characteristics (Hawley & Stone 1995). 
We employed the massively-parallel version of the 
code, ZEUS-MP (Fiedler 1998, Norman 2000), to 
produce a data set at resolution 512^. 

We used the uniform driving mechanism de- 
scribed in Mac Low (1999). At each timestep, a 
fixed pattern of velocity perturbations is added. 



with the amplitude adjusted such that the energy 
input rate is kept constant. The driving results in 
an rms Mach number of = 5 for models of se- 
ries Q and of = 10 for such of scries £ (sec Tabic 
1). The mechanism for generating the perturba- 
tion field allows us to select fixed spatial ranges. 
All models presented in this work employ driving 
at wavenumbers k=l — 2 waves per box length L. 

Our measurements begin at system time t = 
0.0, when the model has reached an equilibrium 
state between the energy dissipation rate due to 
shock interaction and numerical diffusion, and the 
driving energy input rate. Self-gravity is imple- 
mented via an FFT-Poisson solver for Cartesian 
coordinates (Burkert & Bodenheimer 1993). It is 
activated as soon as the model reaches the dissi- 
pation equilibrium at t = 0.0. 

The isothermal equation of state renders the 
system scale free. In code units, the length of the 
box is L = 2 and the mass is M = 1. We dimen- 
sionalize the results by choosing the Jeans length 
Aj, Jeans mass Mj, and free fall time tff (see 
Klessen, Heitsch, & Mac Low 2000 for discussion 
of the scaling). The number of Jeans masses and 
Jeans lengths in each run are given in Table 1. 

All the models arc initialized with a uniform 
magnetic field stretching across the box along the 
2:-direction. The field becomes distorted over time, 
but the magnetic flux $ through the boundaries 
should be constant according to the equations of 
ideal MHD, and is typically preserved by the code 
to a relative accuracy of 10~^ - lO""*^. 

As in HMK, we scale the initial magnetic field- 
strength by the plasma beta f3 = SttP/B"^. The 
parameter /3 is directly related to the critical mag- 
netic flux required to prevent global gravitational 
collapse (see McKee et al 1993) 

where the constant c$ is found to be 0.13 for a uni- 
formly magnetized sphere (Mouschovias & Spitzer 
1976) and 0.16 for a uniformly magnetized sheet 
(Nakano & Nakamura 1978). Here we choose the 
latter value. We find 




where M/Mc = $c/^ is the ratio of mass to crit- 
ical mass as in HMK. The models discussed in 
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this paper span the range (3 — 0.05 — 4.04 and 
M/Mc = 0.88 - 8.3. The Alfven Mach number 
A^A is related to (3 and the sonic Mach number 
M by Ma = M{(3/2y/^. Thus, we consider both 
sub-Alfvenic and super- Alfvenic models. 

A detailed discussion of the model sets is given 
by HMK. Here we describe their physical proper- 
ties just briefly. The supersonic turbulence sup- 
ports the gravitationally unstable region against 
global collapse, however, it cannot prevent local 
collapse, even in the presence of magnetic fields 
too weak to provide magnetostatic support. Col- 
lapsing regions evolve from shock-induced fila- 
ments and cores. The models reach density con- 
trasts of 2— 3 orders of magnitude above and below 
the mean density. The gravitationally bound re- 
gions can be interpreted as the initial stages of a 
protostellar core which may subsequently evolve 
into stars. Once a core begins to collapse, we can- 
not resolve it well enough to follow its evolution 
further (Truelove et al. 1997, HMK). 

The isotropy of a model depends largely on its 
field strength. For the weak-field-model £hla, the 
resulting distribution of magnetic energy is fairly 
isotropic, whereas for the models with stronger 
fields (e.g. £hld), the field imprints its initial di- 
rection onto the gas flows. The flows, in turn, act 
on the field. Two quantities of interest for the po- 
larization studies are the ratio of the mean field 
energy to the total magnetic energy, (5)^/(5^), 
and the dispersion in the angle S between the local 
magnetic field direction and the mean direction, 
cr^((5) (the two are not equivalent; for example, 
the energy in the field could be increased by col- 
lecting it into strong unidirectional filaments with 
no dispersion in angle). Both quantities are given 
in Table 2. 

Equipartition between turbulent kinetic and 
turbulent magnetic energy is often assumed in 
astrophysics; it has been shown to hold rigor- 
ously only in certain cases, such as weak Alfven 
wave turbulence (Zweibel & McKee 1995) and the 
incompressible Alfvenic turbulence modelled by 
Goldreich & Sridhar (1997). Observations do not 
yet yield clear evidence for or against equiparti- 
tion, although data collected by Crutcher (1999) 
tend to speak against it. The models discussed 
here are not in exact equipartition; the ratio 
^ = E^^^/Ekin of turbulent magnetic to turbu- 
lent kinetic energy for each model is listed in Ta- 



ble 2, and is typically a few tenths. It is unclear 
whether these departures from equipartition oc- 
cur for physical reasons related to the nature of 
nonlinear, compressible MHD turbulence or occur 
because of numerical diffusivity or the nature of 
the forcing. We take the deviation from equipar- 
tition into account in § 3.3 to correctly interpret 
our polarization maps. 

2.2. Generating Polarization Maps 

The model cubes include density, magnetic 
fields and velocity fields. These data allow us 
to derive the Stokes parameters Q and U, and the 
polarized intensity P according to Zweibel (1996) 

P = Q + ^U^Jfiy)^-^0^dy (3) 

where and give the field vectors in the plane 
of sky perpendicular to the line of sight and are 
taken directly from the simulations. We integrate 
along the line of sight in y-direction. 

The function f{y) is a weighting function which 
accounts for the density, emissivity, and polarizing 
properties of the dust grains. In this paper, we 
take it to be the gas density normalized by the 
mean density. The polarized intensity is \P\ = 
\JQ'^ -\- U^, and the polarization angle is 

(/>= iarctan^. (4) 

Equation (3) is an approximate solution of the 
full radiative transfer equation for the Stokes pa- 
rameters (Martin 1974, Lee & Draine 1985), valid 
for small polarization and low optical cleptli. At 
far infrared wavelengths, and at typical column 
densities for molecular clouds, the medium can 
safely be assumed to be optically thin (Hildebrand 
ct al. 2000). 

Polarization maps are made at finite resolution. 
We simulate the telescope beam by applying a 
Gaussian filter of the form 



F{x,x') = 




to the complex polarization P given in equation 
(3). The width w of the smoothing filter should 
not exceed one eighth of the box length L. For w > 
L/8, the tails of the filter - exceeding the map's 
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area would contribute to such an extent that 
neglecting them would yield a too small average. 
For each model, we generated polarization maps 
with a set of smoothing widths < w < L/8. 

3. Results 

The 2D polarization maps serve a threefold pur- 
pose: We discuss how their structure depends on 
self-gravity and limited observational resolution 
(§ 3.1), and we show that there is no preferred 
alignment between filamentary structures in our 
simulations and the magnetic fields (§3.2). Fi- 
nally, in § 3.3 we discuss the Chandrasekhar- Fermi 
method of determining the mean field and an ex- 
tension, which determines the rms field, both of 
which we test and calibrate. 

3.1. Structure in Polarization Maps 

The polarization maps showing column den- 
sity and polarization vectors (Fig. 1 left column, 
model £hld) are highly structured. Shock fronts 
moving through the gas initiate formation of fila- 
ments and knots. After one free-fall time (lower 
left panel in Fig. 1), the filaments fragment and 
concentrate due to self-gravity. Qualitatively, no 
influence of self-gravity on the large-scale struc- 
ture of the field is discernible, although there is 
some effect on the smallest scales. 

Smoothing these maps with a Gaussian filter of 
4 pixels width (Fig. 1 right column, £hld) loads 
to a clumpier appearance of the previously well- 
defined filaments. Any substructure in these is 
lost. Single shock structures are smeared out. As 
expected, the field appears more uniform, thus in 
our case indicating more and more its initial orien- 
tation. Small scale variations in the field indicat- 
ing turbulence and defining the turbulent cascade 
are lost due to the smoothing. 

Figure 2 quantifies this effect in model Shld. 
It shows the power spectrum of the angle between 
the local and the mean magnetic field. This is 
equivalent to the line of sight averaged spectrum 
of the magnetic field fiuctuation amplitude (recall 
that the total power in fluctuations is given in Ta- 
ble 2). Increasing the smoothing width from w = 
to w = L/8 results in a power loss of 96%. Struc- 
tures at a wave number of fc = 8 are suppressed 
by more than two orders of magnitude. The dia- 
monds in Figure 2 indicate the power spectrum of 



perturbed against mean field energy E^^g/Emag, 
where E^^^ is the magnetic energy corresponding 
to the perturbed field components perpendicular 
to the mean field. We use this as a measure of the 
true disorder in the magnetic field and as a gross 
check whether the polarization angles mirror the 
behaviour of the energies. Some power loss be- 
tween the true angular dispersion and even the 
unsmoothed fluctuation amplitude is inevitable, 
because the latter involves line of sight averaging 
while the former is a true 3D quantity. 

Inspection of Figure (1) shows no dramatic ef- 
fects caused by self gravity after one free-fall time. 
Figure (3) shows that the width of the distribu- 
tion of polarization angles does not change signifi- 
cantly under the effect of self-gravitation, but does 
certainly decrease with increasing beam width, as 
expected. The emerging asymmetry results from 
the reduced statistics due to the increasing beam 
width. Although every pixel contributes to the 
histogram, the number of independent measure- 
ments decreases with increasing beam width. 

3.2. Shock-induced Filaments 

There is an ongoing debate on the mechanisms 
generating the observed filamentary structures of 
molecular clouds (e.g. Loren 1989, Johnstone & 
Bally 1999, Matthews & Wilson 2000). Shocks 
come to mind as a natural explanation, either 
compressing the material (as in the simulations 
presented) or generating downstream flows result- 
ing in the filaments (Loren 1989). The alignment 
of the magnetic field with the filaments has been 
used as a test for the validity of different fila- 
ment models. The observations seem to favour no 
definite alignment. Maps of OMC-III (Matthews 
& Wilson 2000) show a perpendicular alignment, 
thus perhaps supporting the model of Fiege & Pu- 
dritz (2000a, 2000b), in which they propose helical 
fields to confine filamentary structures. Rizzo et 
al. (1998) found parallel and perpendicular align- 
ment in their maps of background starlight polar- 
ization for Lupus 1 and Lupus 4, whereas Rao et 
al. (1998) find varying alignments for Orion-KL. 
Further examples of varying field alignments can 
be found in Dotson (1996) and Ward- Thompson 
et al. (2000). 

In the simiilations presented here, the filaments 
are solely due to shock interactions. We do not 
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find a preferred alignment of the magnetic field 
with the filaments. We illustrate this in Figure 4, 
which shows the full polarization map for the 512^ 
zone model Srla, and with a selection of filaments 
of the same model in Figure 5. A detailed study 
of the filaments will be presented in a forthcoming 
paper. 

3.3. Estimating the Field Strength: Chan- 
drEisekhar- Fermi-Method 

Chandrasekhar & Fermi (1953; CF) suggested 
a method of estimating the mean magnetic field 
strength (B) in the galactic spiral arms. It relates 
the line of sight velocity dispersion cr(w;os) to the 
dispersion of polarization angles a{taii6) around 
a mean field component. The CF method rests 
on two main assumptions: that the fluctuations 
are isotropic about the direction of the mean field, 
and that there is equipartition between the mean 
kinetic and magnetic energies in the fluctuations. 
Under these conditions, the mean field component 
{B) is given by 



{Bf = 477/9 



(6) 



cr(tan5)^ 

where p stands for the mean density. Because 
of the angular variations in the denominator, the 
result depends strongly on the actual mean field 
strength, and will in fact only be meaningful if 
there is a noticable mean field component. In 
terms of the parameter ^, the ratio of turbulent 
magnetic to turbulent kinetic energy listed in Ta- 
ble 2, equation (6) generalizes to 



iBf = 4.p4^e 



(7) 



(7 (tan (5)^ 

We initially set ^ = 1 when implementing the 
CF-method, because observations do not yield 
clear evidence for or against equipartition (e.g. 
Crutchcr 1999), but we correct for ^ ^ 1 later 
(see Tabic 2). 

In the following paragraphs, we present a modi- 
fication of the CF-mcthod (§ 3.3.1), compare both 
methods with the help of model data (§ 3.3.2) and 
discuss the effect of limited telescope resolution on 
the resulting field strength estimates in § 3.3.4. 

3.3.1. A Modification 

The CF-method in its original form estimates 
the mean magnetic field. When the mean field 



is much less than the rms field, the dispersion in 
fluctuation angle a^{tand) is dominated by points 
where 6 « 7r/2. Since {B) oc a{ta.n6)~^ (see 
eq. (6)), the result is usually an underestimate of 
(B). 

The CF-method can be extended to yield an 
estimate {B'^) which is free of this problem. Sup- 
pose {B) = z{B) and the x — z plane is the plane 
of the sky. Then 



'B) 



{l^Bl) + {Bl) + {Bl). (8) 



By definition of the angle 5 

B^ = (B) tan 5. 



(9) 



Squaring and averaging equation (9), assuming 
that the magnetic field fluctuations energies are 
the same for all components (see § 3.3.6), and us- 
ing equation (6), equation (8) becomes 

^^'^ = ^"""^^^W^^ + 3<7(tan5)^). (10) 

When the fleld is highly disordered, and the 
dispersion in polarization angles is large, equation 
(10) reflects the underlying physical assumptions 
of equipartition between kinetic and magnetic en- 
ergy together with isotropy. The formula is valid 
for arbitrary ratios of turbulent magnetic to tur- 
bulent kinetic energy if we multiply the RHS by 
^. When the field is highly ordered, comparison 
of equation (10) with equation (6) shows that the 
mean and rms fields are nearly the same. 

3.3.2. Calibration with Model Data 

We plotted the ratio of the estimated and model 
field strength a = B^^st/ B^nodei in Figure 6. The 
upper left panel displays a for the CF-method, 
averaged over two lines of sight perpendicular to 
the mean field direction and over three physical 
times. The lower left panel shows the dispersion 
of a single measurement relative to the correspond- 
ing mean value. It decreases with increasing field 
strength. For strong fields, the CF-method over- 
estimates the field strength derived from the un- 
smoothed maps (star symbols in Fig. 6) by a fac- 
tor of 2 to 3. Thus we confirm the claim made by 
Ostriker et al. (2001). For weaker fields however, 
the method starts to develop a significant scatter, 
as shown in Figure 6 (lower panel). 
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The right cohimn in Figure 6 show the corre- 
sponding results for the CF-method in its modified 
version according to equation (10). For the un- 
smoothcd case (again star symbols), the strongest 
field is overestimated by a ~ 3. Whereas the orig- 
inal method breaks down at a field strength of 
Bmodei = 1-2, the modified version still yields ac- 
ceptable results for Bmodei = 0.9. 

The minimum field strength alone however 
does not tell us much about the reliability of the 
method. The parameter of interest is the Alfven 
Mach number A4a = (w^)^^^/c^, as the angular 
variations a{ta,n6) of polarization not only depend 
on the energy content in the field, but on the tur- 
bulent kinetic energy as well. A field strength of 
Bmodei = 0.9 would correspond to Ma = 1-25, in 
the parameter set of models Q. With typical pa- 
rameters of Cg = 200m/s and n(H2) = 10^ cm~^, 
we can then scale Bmodei = 7nG. Field strengths 
in molecular clouds are seen to be mostly larger 
than B > lOfiG (Crutcher 1999). In denser re- 
gions, the field strengths seem to increase with 
density as JB oc (Crutcher 1999), whereas 

the nonthermal line widths (the "turbulence" ) de- 
creases with decreasing size as Avnt oc 
(Caselli & Myers 1995). Thus we conclude, that 
both methods should yield acceptably accurate 
results for molecular cloud regions up to their 
densest parts, as long as the observations sample 
the angular variations sufficiently (see § 3.3.4). 

We have to address the question of energy 
equipartition. Since ^ is generally less than unity 
in our models (see Table 2), the CF-method and 
its extension overestimate both (B) and (B^) by 
a factor of l/V^. Figure 7 shows a = Besi./ Bmodei 
"corrected" for non-equipartition. Both methods 
now hit the model field strength at a factor be- 
tween 1 and 1.5, at least for sufficiently large field 
strengths and for unsmoothed maps. The correc- 
tion slightly reduces the largest deviations of the 
weak-field estimates. We conclude, that for a suffi- 
ciently well resolved (sec below) polarization map 
and for large field strength, both methods yield 
reliable results. A main uncertainty factor is the 
ratio of magnetic to kinetic energy ^ in the ob- 
served region. 

Whereas the original CF-method underesti- 
mates low field strengths, the modified version 
overestimates them. For weak fields, the polariza- 
tion angles can reach 90° with respect to the mean 



field, which yields |tan^| = oo in equations (6) 
and (10). 

As the original CF-method, the modified ver- 
sion yields overestimated field strengths with in- 
creasing smoothing beam width. From the lower 
row in Figure 6 we conclude that even for the 
strongest fields the relative scatter is of order 20%. 
Smoothing the maps leads to smaller scatter ex- 
cept for the weakest field strengths. 

The effect of self-gravity on the field strength 
estimates is minute (Fig. 8). For small smoothing 
widths w, the varying small scale structure due to 
self-gravity leads to some statistical scatter, which 
is averaged out for large w. The CF-methods seem 
to be insensitive to effects of self-gravity on small 
scales. 

3.3.3. A Recipe 

The fact that acF and Urms vary by roughly the 
same logarithmic range, but with opposite signs, 
suggests that by taking the geometric mean 

= 4.p4^(l+3.(tan^)^)V^ ^''^ 

(T(tanc))^ 

we would arrive at a more accurate estimate for 
the actual field strength. Applying this recipe 
yields Figure 9 (corresponding to Fig. 7). Clearly, 
the large deviations of acF and arms cancel suf- 
ficiently to yield estimates accurate to a factor 
of Ogm ~ 2.5 even for the weakest, unsmoothed 
fields, and is more accurate than that for moder- 
ately strong fields. Smoothed maps of course lead 
again to a systematic ovcrcstimation, here up to 
a factor of 10. We would like to emphasize that 
equation (11) is only motivated by the (logarith- 
mically) comparable deviations of acF and arms- 
To test this recipe, we applied it to a set of 
models with varying physical properties (Fig. 10). 
Models LI, L2 and L3 are a time scries of decaying 
turbulence taken from Mac Low (1999). Models 
MC81, MA81 and MC41 (Mac Low et al. 1998) 
are simulations of driven turbulence as series Q, 
but with driving wave length of fc = 3 — 4 (MC41) 
and A; = 7 - 8 (MA81, MC81). See Tables 1 
and 2 for their parameters. All these models have 
weak initial fields, corresponding to Qilb. As in 
Figure 9, the scatter of agm is more than one order 
of magnitude lower than for acF and arms- The 
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unsmoothcd estimates are accurate again up to 
a factor of agm ~ 2.5. Again, smoothed fields 
are likely to be overestimated. Note that 

^rms on 

the whole leads already to slightly more accurate 
estimates than acF- 

We conclude that the geometric mean Ugm gives 
more reliable estimates than acF or agm even for 
very weak fields, being more accurate as well for 
moderate and strong fields. Thus, this recipe may 
prove useful when being applied to observational 
data. 

3.3.4- Smoothing, subsampling and their effect in 
a turbulent medium 

Field estimates derived from smoothed maps 
tend to overestimate the field strength with in- 
creasing "beam width" (Fig. 6, all symbols except 
stars). We will discuss how limited telescope res- 
olution in polarization maps of turbulent regions 
can lead to this overestimation. 

Both CF-mcthods rely on perturbations of the 
magnetic field around a mean component. If all 
the power of these perturbations resided on spa- 
tial scales 1/k = w < w, increasing the smoothing 
width above w would have no effect on the power 
distribution of polarization angles, and varying 
smoothing widths would yield identical results. 
On the other hand, if the perturbations grew 
with increasing length scale, as e.g. in a turbu- 
lent power spectrum, a larger smoothing width 
w would result in a greater power loss on larger 
scales, and thus to an overestimation of the field. 
Note that in equations (6) and (10) we find the 
angular variations in the denominator. So losing 
power in the variations increases the field estimate. 
The turbulence is driven between wave numbers 
k = 1 — 2 in our simulations, thus letting the code 
evolve a turbulent cascade down to the dissipation 
scale at grid cell size. This turbulent cascade in 
our models is well reflected in the power spectrum 
of angular variations a{d) (Fig. 2). We lose more 
and more power when applying larger smoothing 
widths w. This power loss should depend on the 
steepness of the power spectrum. This suggests 
that it might be possible to measure the spectrum 
of magnetic field fluctuations by smoothing polar- 
ization maps by different filter widths, and com- 
paring the dispersion of polarization angles. 

To illustrate this relation, we will describe the 



smoothing as a convolution and then will deter- 
mine analytically the loss of power in a spectrum 
due to smoothing. The convolution would be 

/oc 
A{x)F{x - x') dx (12) 
-oo 

with the wave function A{x) 

A{x) = 51 ( r ) sin(27rfca;) (13) 

fe=i ^ ^ 

and the Gaussian filter function F{x — x') as in 
equation (5). We identify A{x) with the Bl/{B)'^ 
of equation (9) . Instead of convolving directly, we 
use the Fourier transforms J^[A\ and J^[F\ of equa- 
tions (13) and (5) to derive the power spectrum of 
the convolved function via the convolution 

theorem 

A = AF, (14) 

where the hatted quantities denote Fourier trans- 
formations. The Fourier transform of a Gaussian 
yields again a Gaussian, in our case with an addi- 
tional phase shift. The frequency space variable is 

U). 




(15) 



Applying the Fourier transform to the wave func- 
tion A{x) and integrating over k = [1, oo], we get 




Now we can multiply A and F, and taking the 
absolute value gives us the power spectrum of the 
convolved function A. 

|i|=7rexp(-i(a;<) (17) 

As we are interested in the loss of power caused 
by the smoothing, we now try to determine the 
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amount of power remaining in the spectrum after 
smoothing, i.e. the quantity 

K{w) = J TT cxp ^—^{ujw)^^ ^ 

(18) 

For 7 < 1, lim^^Q K{w) = oo. Figure 11 illus- 
trates the behaviour of K{w) for several values of 
7. Note that this is a log-linear plot, so that, for 
example, K{w,j = 3/4) cx log(w)^. Curves for 
7 > 1 have negative, the ones for 7 < 1 positive 
curvature. Normalizing the curves for 7 < 1 with 
K{0) and plotting them together with the inverse 
of the ratio of estimated to model field strength 
acF of model Ghlh yields Figure 12. We normal- 
ized acF to w = as well, in order to compare 
the measurements to the analytic curves K{w). 
We note that the measurements start to flatten 
off for small w/L because of line of sight aver- 
aging, which we do not include in the analytic 
derivation. The same argument applies for ob- 
servations. Thus, even in the ideal case of w = 0, 
the measurement may not reach a = 1 but may 
overestimate the actual field strength. Of course, 
this would crucially depend on the length of line 
of sight and, generally, on optical depth effects. 
They, however, are neglegible in the far-infrared 
regime. Because we need a certain length in the 
line of sight in order to derive a meaningful ve- 
locity dispersion a{vios), this overestimate due to 
line-of-sight averaging gives an intrinsic error in 
the CF-methods. The power spectrum of angu- 
lar variations in the simulations is by no means a 
clear power law, due to numerical dissipation and 
the signatures of the driving scale. Thus, the mea- 
sured l/acF (Fig- 12) does not correspond clearly 
to one of the analytic curves. However, if observa- 
tions show a well defined power law in the angular 
variations a (tan 6), and provide sufficient resolu- 
tion to create a series of maps smoothed over larger 
and larger areas, then it should be possible to ex- 
tract information about the exponent of the power 
law and to correct the field estimate as well. Once 
a power law description is found for the observed 
region, K{'w)/K{0) supplies a correction factor for 
the field strength estimate B^st- 

In polarization maps derived from observations, 
an additional effect could potentially lead to a sys- 
tematic overestimation. When we derive the field 
estimates from our models, we use the complete 



domain, thus including all wave lengths occurring 
in the problem. However, it is not clear whether 
the observations do trace the largest scales in the 
field perturbations. If the power of angular vari- 
ations is distributed in a power spectrum which 
decreases with increasing wave number, we will 
lose overproportionally much power when estimat- 
ing the field strength within a domain which is 
smaller than the largest wave length contribut- 
ing to the perturbations. Figure 13 illustrates 
this effect contributing to a systematic overesti- 
mation of the field strength. It shows the ra- 
tios a ~ B^st/ Bmodei for subframes of 32^, 64^ 
and 128^ cells for model Ghlh at 256^ resolution. 
We averaged the estimates from all possible sub- 
frames. The error bars indicate the standard devi- 
ations of the means. With decreasing frame size, 
the field is overestimated slightly, thus indicating 
a decrease of a{tan5). However, bearing in mind 
the large scatter around the mean values in Figure 
13, this eff'ect is not very distinct. 

We conclude that, although the CF-methods 
are applicable in principle to molecular clouds, 
the resulting field strengths are most likely to be 
systematically overestimated for two reasons: (a) 
Any smoothing introduced by the finite telescope 
beam size will lead to a systematic removal of 
large-scale perturbations and thus to an overes- 
timation of the field strength, (b) Deriving field 
estimates for a domain smaller than the maximum 
perturbation wave length present leads to power 
loss on large scales and thus as well to an overes- 
timated field strength. This systematic overesti- 
mation due to limited resolution is different from 
uncertainties because of non-equipartition as dis- 
cussed in § 3.3.2. In order to draw these conclu- 
sions we assume that the imderlying power spec- 
trum of perturbations follows a power-law-like dis- 
tribution, at least, that it decreases with increas- 
ing wave number. 

3.3.5. Resolution Test 

Figure 14 presents a test for resolution with 
models Qilh and Qhlh at 128'' and 256^ grid cells. 
The beam width W256 is given in pixels. Note that 
in order to compare both models, we have to com- 
pare them at the same physical beam width, i.e. 
we have to shift a data point of model GHh at 
Wi28 = 4 such that it corresponds to a point of 
W256 = 8. We corrected all estimates with 1/^ (see 
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Tabic 2) in order to account for non-cquipartition. 
Discrepancies at small scales are to be expected 
as model ^il/i then already reaches the grid scale. 
Figure 14 allows us to conclude that numerical 
resolution does not affect our calibration of the 
CF-methods. Thus we feel justified to draw con- 
clusions for observations of turbulent regions from 
our simulated observations of computational tur- 
bulence. 

3.3.6. Isotropy of simulated turbulence 

We still have to substantiate the claim leading 
to equation (8), namely, that our simulated tur- 
bulent flows are isotropic enough to validate that 
the field components perpendicular to the mean 
field (B) = z{B) obey {B^) = {By) and that 
the field perturbations in all three directions ful- 
fiU (Asf) = {Bl) = {Bl). The upper panel of 
Figure 15 supports the second assumption. There 
occur only minor differences in the rms values. 
Note that the initially homogeneous background 
field is oriented along the ^-direction. Thus, we 
had to subtract this background in order to de- 
termine the corresponding mean values, yielding 
{Bz) (lower panel of Fig. 15). The scatter of the 
means is considerably larger, however their abso- 
lute values are three to four orders of magnitudes 
lower than the corresponding field strengths. Ide- 
ally, the means should vanish, especially in the di- 
rections perpendicular to the initial field, namely 
the X- and y-direction. Deviations from this ideal 
case are probably due to numerical diffusion in the 
code. This results in not conserving the total flux 
per coordinate direction. The small shifts in the 
means thus indicate how well the code actually 
does conserve the total flux. 

4. Conclusions cind Implications 

Reliable methods for measuring and mapping 
magnetic fields in molecular clouds are urgently 
needed in order to assess the role of magnetic 
fields in the dynamics of molecular clouds and in 
star formation. Probes based on far infrared po- 
larimetry of the thermal emission from magneti- 
cally aligned dust are now in widespread use. Far- 
IR polarimetry is complementary to the Zeeman 
method, which is the traditional means of measur- 
ing the field. There arc few independent checks of 
its accuracy as a diagnostic of magnetic fields. 



In this paper we used numerical simulations 
of turbulent, magnetized, self gravitating molec- 
ular clouds to assess the reliability of polarization 
maps. Although such simulations qualitatively re- 
produce many of the features seen in molecular 
clouds, we do not expect the computed densities, 
velocities, and magnetic fields to be completely re- 
alistic. We undertook the project with the expec- 
tation that analysis of simulated data would bring 
out the same issues as analysis of real data. 

We assumed that the polarized emission is op- 
tically thin and proportional to gas density. The 
first assumption should be excellent. The second 
could be modified, for example by modelling the 
grain temperature point to point within the cloud. 
It might be possible to pick out particular struc- 
tures along the line of sight by comparing maps 
at different wavelengths, but establishing this is 
beyond the scope of the present paper. 

The salient properties of the models are given 
in Tables 1 and 2. Based on these models, our 
principal conclusions are as follows: 

1. Filaments produced by shocks do not show a 
preferred alignment with the magnetic field. 

2. Self- gravity has no discernible effect on the 
structure of the magnetic fields in our simu- 
lations. However, we have to bear mind that 
our self-gravitating regions are small with re- 
spect to the simulated region. 

3. The Ghandrasekhar-Fermi method in its 
original version yields field estimates ac- 
curate up to ocF = Bcf/{B) = 2. a 
modified version determining the rms field 
B^f results in in a™« = B^f/Brm. = 2. 
This version proves to work more reliably 
for slightly weaker field strengths. Both 
methods should be applicable to magnetic 
fieldstrengths typical of molecular clouds. 

4. The geometrical mean agm = y/acF 
yields field estimates acciirate up to a fac- 
tor of w 2.5 even for the weakest fields. 
For moderate and stronger fields, agm ~ 
1. Thus, this recipe improves the reliabil- 
ity of field strength estimates for the whole 
range of physically realistic fields in molec- 
ular clouds, and may prove useful when ap- 
plied to observational data. 
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5. The original and modified Chandrasekhar- 
Fermi methods are based on the assump- 
tion of equipartition between turbulent ki- 
netic and magnetic energy. The field fluctu- 
ations in the models are below equipartition, 
leading to an overestimate of the field. 

6. Limited angular resolution again leads to 
systematic overestimation of the magnetic 
field strength for both methods. Two effects 
can cause this overestimation: (a) In a power 
spectrum decreasing with increasing wave 
number, smoothing leads to an overpropor- 
tionally large power loss on larger scales. 
Thus the angular variations are underesti- 
mated, which leads to larger field strengths, 
(b) If the domain used for deriving the field 
estimate does not trace the largest wave 
modes, this again leads to power loss on the 
largest scales. 

7. If the power spectrum of magnetic field fiuc- 
tuations is a power law, field estimates made 
from a series of maps with different smooth- 
ing widths can in principle be used to esti- 
mate the exponent of the fluctuation spec- 
trum, assuming that this is well enough de- 
fined over at least two decades. With this 
information, one could correct for the over- 
estimated field strengths. 
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Table 1 
Model list 



Name 


Resolution 


^drv 


P 


M/Mer 


Aj 


nj 


Srla 


5123 


1 - 2 


4.04 


8.3 


0.501 


64 


£hla 


2563 


1 - 2 


4.04 


8.3 


0.501 


64 


Shld 


2563 


1-2 


0.20 


1.8 


0.501 


64 


Qila 


1283 


1 - 2 


4.04 


8.3 


1.067 


6.5 


Gilb 


1283 


1 - 2 


1.13 


4.4 


1.067 


6.5 


Qilc 


1283 


1 - 2 


0.50 


2.9 


1.067 


6.5 


gud 


1283 


1 - 2 


0.18 


1.7 


1.067 


6.5 


gne 


1283 


1 - 2 


0.13 


1.5 


1.067 


6.5 


Qilf 


1283 


1 - 2 


0.09 


1.3 


1.067 


6.5 


gug 


1283 


1 - 2 


0.07 


1.1 


1.067 


6.5 


gnh 


1283 


1 - 2 


0.05 


0.8 


1.067 


6.5 


ghlh 


2563 


1 - 2 


0.05 


0.8 


1.067 


6.5 


L 


2563 


1 - 8 


2.00 


00 


(...) 


(...) 


MCSl 


1283 


7-8 


2.00 


00 


(...) 


(...) 


MAil 


1283 


7-8 


2.00 


oo 


(...) 


(...) 


MCAl 


1283 


3-4 


2.00 


oo 


(...) 


(...) 



Note. — /? = Pth/Pmag = 8'jTc'^p/B'^. M/Mcr gives the ratio of cloud mass to critical mass according to 
equation (1). Aj is the Jeans length, with a total box side length oi L — 2, and nj is the number of Jeans 
masses in the box. The non-selfgravitating models L, MC81, MASl and MCAl are taken from Mac Low et 
al. (1998) and Mac Low (1999). 
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Table 2 

Characteristic parameters of all models 



Name 


(3 




Efnag 1 ^kin 




^{5) ["] 


Srla 


4.04 


0.96 


1.03 


13.06 


22.3 


Shla 


4.04 


0.35 


0.37 


12.00 


15.0 


Shld 


0.20 


0.52 


1.30 


0.79 


14.0 


Qila 


4.04 


0.21 


0.23 


8.90 


20.0 


gub 


1.13 


0.36 


0.45 


3.90 


17.0 


guc 


0.50 


0.69 


0.91 


3.10 


13.0 


gild 


0.18 


0.48 


1.00 


0.88 


14.0 


gne 


0.13 


0.31 


1.00 


0.45 


13.0 


guf 


0.09 


0.28 


1.20 


0.32 


11.0 


gug 


0.07 


0.27 


1.30 


0.27 


11.0 


gnh 


0.05 


0.24 


1.70 


0.17 


8.7 


ghlh 


0.05 


0.26 


2.00 


0.15 


9.1 


LI 


2.00 


0.83 


1.15 


2.55 


15.4 


L2 


2.00 


1.07 


1.78 


1.53 


14.4 


L3 


2.00 


1.03 


2.45 


0.73 


12.7 


MC81 


2.00 


0.41 


0.50 


4.45 


16.4 


MA81 


2.00 


0.79 


1.31 


1.51 


14.5 


MCAl 


2.00 


0.39 


0.48 


7.82 


15.4 



Note. — All values are taken at t = 0.0 except for models LI to i3, who are a time series of decaying 
turbulence (Mac Low 1999). C = ElZg/Ektn, the ratio of turbulent magnetic to turbulent kinetic energy, 
where E^^g does not include the mean (z-)field contribution. Emag/Ekin is the ratio of total magnetic 
energy (including uniform component in 2:-direction) to kinetic energy, Bf^^^/{B)^ is the ratio between the 
turbulent magnetic energy (= E^^^g) and mean field energy, and a{5) gives the dispersion of field angle 
around the mean field direction in degrees. 
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Fig. 1. — Le/t column: Surface density (corresponding to continuum density) of the full computational 
domain, overplottcd with polarization vectors as determined by equation (3) for the sclf-gravitating model 
£hld. The magnetic field is supercritical with M/Mc = 1.8. At t = 0.0 (upper panel), gravity has been 
switched on, which shows its effect in the lower panel, at t = 1.0, a free-fall time later. The initially uniform 
field is vertically orientated (as well as for all other maps shown). Right column: Surface density and 
polarization vectors as in left column, but smoothed with a Gaussian filter of w = 4 pixel. The field gets 
more ordered and the previously well discernible shock structures are blurred or even come out as clumpy, 
filamentary structures. 
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Fig. 2. — Power spectrum of polarization angles with respect to mean magnetic field of model Shld and 
ratio AEmag/ Emag, Corresponding to {B"^) / (By)'^ , the perturbed against mean field (diamonds) against 
wavenumber k. The energy ratio is shifted horizontally by A\og{AEmag / Emag) = 0-5 for clarity, w denotes 
the smoothing width, < w < L/8, with L = 256 the box length. At w = L/8, 96% of the power in 
polarization angles is lost. 
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Fig. 3. — Distribution of polarization angles for two times and various smoothing beam widths w of model 
£hld. The bottom panel starts with t = 0.0, when gravity is switched on. The upper panel shows the angle 
distribution one free-fall time later. With increasing w, the distributions get narrower, whereas self-gravity 
does not affect the distributions significantly except for the smallest scales {w = Q). 
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Fig. 4. — Surface density of the full computational domain of 512^ cells, overplotted with magnetic polar- 
ization vectors as determined by equation (3) for the not yet self-gravitating model £rla at t = 0.0. The 
magnetic field is supercritical by a factor of 8.3. The white boxes mark the filament locations as shown in 
Figure 5. Every eighth polarization vector is shown. 
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Fig. 5. — Magnetic polarization vectors and surface density for six selected filaments of model £rla as 
denoted in Figure 4. ^0 
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Fig. 6. — Overestimation factors acF = Bcf/{B) and a^ms = B^p"^ / Brms for the CF-method in its original 
and modified version for all models of type Q, averaged over two lines of sight and three physical times (upper 
panels). The dotted lines correspond to a = 1, 2, 3. Symbols stand for smoothing beams applied. Increasing 
beam sizes lead to systematic overestimation of the field strengths. The beam size is given in zones, with a 
total grid size of 128 zones. The lower panels contain the relative standard deviations with respect to the 
corresponding mean value. 
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Fig. 7. — Overestimation factors acF = Bcf/{B) and arms = B^p'^/Brms for the CF-method in its original 
and modified version for all models of type Q as in Figure 6, bnt after correcting for non-equipartition (see 
Table 2) with . The dotted lines correspond to a = 1, 2, 3. Both methods return the model field strength 
nearly exactly for sufficiently high field strength and no smoothing. 
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Fig. 8. Ovcrcstimation factors acF = Bcf/{B), arms = B^p^/Brms and Ugm = {o.cpO'rmsY^'^ for the 
self-gravitating model £hld. At t = 0.0, gravity is switched on. Note that all three panels are plotted in the 
same logarithmic scale. 
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Fig. 9. — Geometric mean = ^/o-cf o-rms against model field strength, averaged and corrected for non- 
equipartition as in Figure 7. (upper panel). Where the separate estimates acF and arms (Fig. 7) miss the 
weak field strengths by a factor > 150 for unsmoothed {w = 0) maps, Ugm yields an estimate accurate up 
to a factor of « 2.5 for the same w = 0-maps. Smoothed fields are still systematically overestimated. The 
lower panel shows the corresponding errors. 
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Fig. 10. — acF, drms and geometric mean agm = \/acF arms for an extended model set (see Table 1, lower 

part), corrected for non-cquipartition as in Figure 7. Where the separate estimates acF and arms scatter 
around a = 1 by up to a factor of 700 for unsmoothed {w = 0) values, Ugm reduces this scatter to a factor 
< 10. 
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Fig. 11. Analytic results for the amount of power remaining in the spectrum K{w) (sec eq. (18)) for five 
values of the power law exponent 7. In this log-linear plot, exponents of 7 > 1 and 7 < 1 can be distinguished 
by their curvature. L gives the box length. 
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Fig. 12. Analytic results for the amount of power remaining in the spectrum K{w) (see eq. (18)) for three 
values of the power law exponent 7. The stars denote the inverse of acF taken from model Qhlh. All curves 
have been normalized to w = in order to compare them. L is the box length. 



27 



cd 



cd 



6 

4 



w = 0.0 



w = 4.0 



<> 







100 200 
framesize N 



100 200 
framesize N 



Fig. 13. — Ovcrcstimation factors acF = Bcf/{B) and a„„s = B'^p'^ / Bmis against linear size of subframes 
for model Qhlh. acF and arms were determined by averaging over all possible subframes for a given box 
size. 
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Fig. 14. — Overestimation factors acF = Bcf/{B) and arms = B^p / Brms against beam width w for the 
strong-field models Qilh and Qhlh at resolutions of 128^ and 256^ cells. The beam width W25q is given in 
pixels, where ^256 refers to the beam width of model Qhlh. We have to compare the different resolutions at 
the same physical beam width, which is why the data points of model GHh are shifted by a factor of 2, i.e. 

Wi2S = 2 ^ W256 = 4. 
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Fig. 15. — rms field values (upper panel) and mean values (lower panel) of the magnetic fields in the three 
coordinate directions x, y and z. The initially uniform field is parallel to the 2;-direction. For the ^-direction, 
we determined the mean value after subtraction of the background field. 
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